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Mathematical regularisation of the nonlinear terms in the Navier-Stokes equations provides a 
systematic approach to deriving subgrid closures for numerical simulations of turbulent flow. By 
construction, these subgrid closures imply existence and uniqueness of strong solutions to the 
^corresponding modelled system of equations. We will consider the large eddy interpretation of 
QQtwo such mathematical regularisation principles, i.e., Leray and LANS— a regularisation. The 
^y-Jjeray principle introduces a smoothed transport velocity as part of the regularised convective 
^^rionlinearity. The LANS— a principle extends the Leray formulation in a natural way in which a 
^1 filtered Kelvin circulation theorem, incorporating the smoothed transport velocity, is explicitly 
^3 " satisfied. These regularisation principles give rise to implied subgrid closures which will be applied 
'■—■in large eddy simulation of turbulent mixing. Comparison with filtered direct numerical simulation 
l/~idata, and with predictions obtained from popular dynamic eddy-viscosity modelling, shows that 
C _ I hese mathematical regularisation models are considerably more accurate, at a lower computational 
***^post. Particularly, the capturing of flow features characteristic of the smaller resolved scales improves 
•-Significantly. Variations in spatial resolution and Reynolds number establish that the Leray model is 
Hmore robust but also slightly less accurate than the LANS— a model. The LANS— a model retains 
C^more of the small-scale variability in the resolved solution. This requires a corresponding increase in 
•the required spatial resolution. When using second order finite volume discretisation, the potential 
^accuracy of the implied LANS— a model is found to be realized by using a grid spacing that is not 
^^latger than the length scale a that appears in the definition of this model. 

Cdr'his paper is associated with the focus-issue Cascade Dynamics: Fundamentals and Modelling. 
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1 Introduction 

1.1 The Large Eddy Simulation (LES) approach for modelling turbulent 
flow 

The need for accurate and efficient prediction of turbulent flow has increased 
in priority during the past few decades, stimulated by the scientific as well as 
the practical significance of this field of interest. Various computational mod- 
elling strategies to address this need have been put forward. Conceptually, 
the simplest approach is to discretise the governing Navier-Stokes equations 
and resolve all dynamically relevant length scales that are contained in an un- 
steady turbulent solution. This direct numerical simulation (DNS) approach 
has proven to be an essential point of departure for understanding fundamental 
properties of turbulence and to provide a reliable under-pinning of theoretical 
and modelling approaches. However, the requirement of numerically resolving 
flow-features down to the Kolmogorov length-scale poses severe restrictions 
on direct numerical simulation. In fact, current computational resources re- 
strict applications of DNS to turbulent flow of modest complexity, i.e., to low 
Reynolds numbers. As a consequence, much research has been directed to de- 
velop simulation strategies that are computationally much less demanding and 
at the same time provide sufficient accuracy for specific applications. This nec- 
essarily implies a coarsening of the flow description which is usually formalised 
in terms of an averaging process that explicitly reduces the dynamic impor- 
tance of the smaller features in a flow, by modelling their effects in terms of the 
resolved features. Different reduced descriptions that have been proposed in 
the literature may be distinguished by their specific underlying averaging pro- 
cess. The well-known Reynolds averaged Navier-Stokes (RaNS) formulation is 
usually defined with reference to either a long-time averaging or an ensemble 
averaging. More recent is the so-called large eddy simulation (LES) which is 
based on the application of a low-pass spatial averaging to the Navier-Stokes 
equations. 

The application of an averaging process to the Navier-Stokes equations im- 
poses an external control over the dynamic importance of the small-scale flow- 
features. Usually, the averaging is supposed to commute with temporal and 
spatial derivatives. This preserves the momentum-conservation form of the 
governing equations, while introducing the so-called 'closure' terms. These 
closure terms are associated with averaging the convective nonlinearity of the 
Eulerian equations and they give rise to momentum fluxes that cannot be fully 
evaluated in terms of the averaged solution alone. In RaNS context, this clo- 
sure problem is defined in terms of the Reynolds stress tensor while in LES the 
turbulent stress tensor arises. In this paper we will restrict ourselves to LES. In 
order to arrive at a meaningful formulation for the prediction of the averaged 
solution, the closure problem must be resolved by introducing a model for the 
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turbulent stress tensor which approximates the dynamic consequences of the 
small-scale flow-features in terms of operations on the averaged solution alone. 
So, while the main virtue of the LES approach is to provide external control 
over the dynamic complexity of the simulation model, its main challenge is to 
provide an accurate closure for the turbulent stresses. 

A significant portion of all LES literature is devoted to developing accurate 
models for the closure problem that are introduced by spatial filtering. Popular 
low-pass filters used in LES, such as the top-hat or Gaussian filters, may be 
characterised by a length-scale parameter that represents the 'width' A of the 
filter. Application of such filters suppresses the contributions of flow- features 
that are more localised than A while the larger, more energy-containing scales 
are virtually unaffected by the filtering. This implies that the closure problem 
for the turbulent stress tensor primarily involves modelling the effects of the 
'sub-A' flow features on the larger flow features. Since in numerical simula- 
tions the filter-width and the spatial resolution are typically chosen to be of 
comparable magnitude, the closure model for the turbulent stresses is usually 
referred to as a subgrid model. 

Various suggestions have been made for obtaining acceptable subgrid mod- 
els. In the absence of a comprehensive statistical theory of turbulence, empiri- 
cal knowledge about small-scale turbulence is essential for the development of 
such subgrid models. This includes a detailed representation of the dissipative 
and dispersive contributions of the sub-A features to the dynamics of the larger 
features. Moreover, guidance for proper subgrid modelling may be obtained 
by requiring the modelled equations to comply as much as possible with ba- 
sic properties of the (filtered) Navier-Stokes equations . In this context one 
may, e.g., consider: preserving Galilean transformation properties of the sub- 
grid model , maintaining realisability conditions in case positive filters are 
adopted [I], accounting for algebraic properties of the turbulent stress tensor 
as basis for dynamic subgrid modelling [3], incorporating resolved small-scale 
information by allowing (approximate) inversion of the LES filter 6 , or com- 
bining filter inversion with the dynamic procedure [7j, just to name a few of 
these 'mathematical' modelling options. 

An extensive account of modelling practices in LES may be found in [H]. 
The variations among the results arising from the different subgrid models 
that have been proposed appears inconsistent with the basic premise of LES 
that small-scale turbulence may be assumed to be 'universal'. Such 'universal- 
ity' would suggest that the characteristic features of the large scales in a flow 
problem should show a degree of insensitivity to the properties of its small-scale 
flow features, provided these features are modelled in the proper universality 
class. Correspondingly, accurate subgrid models should exist which are sim- 
ple, in the sense that one should be able to formulate such models in terms of 
basic properties of the turbulent flow without involving, for example, details 
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related to the geometry of the flow domain. If universality is a proper guiding 
principle for subgrid modelling, then to the extent the results of the numerical 
simulations depend on these modelling steps, the intuitive and empirical mod- 
elling steps that have primarily been considered in the LES literature may not 
belong to the proper universality class. Regardless of whether the premise of 
universality holds, we suggest a alternative systematic approach based on re- 
turning to the turbulence closure problem, concentrating on the mathematical 
and physical aspects of the Navier-Stokes equations, and seeking the simplest 
subgrid model which is still consistent with these fundamental aspects. 

Turbulence simulations based directly on the governing Navier-Stokes equa- 
tions are not only hampered by the complexity of turbulent solutions. In ad- 
dition, there exist outstanding mathematical issues concerned with existence, 
uniqueness and regularity of solutions to these flow-equations. These issues ad- 
dress the fundamental computability of turbulent flow. Thus, they introduce 
uncertainty into the fields of (dedicated) numerical methods and subgrid mod- 
elling of turbulence. In contrast, one may limit possible closure strategies such 
that the final formulation at least is guaranteed to possess a unique solution 
with a priori known smoothness properties. Moreover, it appears natural to 
require that the solutions of the unfiltered Navier-Stokes equations would be 
recovered in the mathematical limit in which the filter-width is sequentially 
reduced to zero. Recently, analysis established such properties for the well- 
known Smagorinsky model [Q|llOj . However, this subgrid model represents the 
dynamics of the small turbulent scales in terms of a nonlinear eddy-viscosity 
only, and this limitation does not do justice to all intricacies of the spatially fil- 
tered convective fluxes. Moreover, from experience gathered in LES using this 
subgrid model, it has become clear that the predictions of the filtered Navier- 
Stokes solution by the Smagorinsky model are in many cases not accurate 
enough to be of practical relevance 

1.2 Averaging, filtering and closure: LES and Lagrangian averaging 
(LA) 

We begin motivating our selection of a class of reduced flow descriptions for 
which existence and regularity properties are available, by recalling a histor- 
ical result for the Navier-Stokes equations which regularised their convective 
nonlinearity. In mathematical analysis, this approach was pioneered in the 
classic work by Leray |12| who introduced a smoothed transport velocity 
in the convective nonlinearity. In detail, Leray replaced u • Vu in the Navier- 
Stokes equations by u • Vu where the over-line on u denotes filtering of the 
velocity u. Both the velocities u and u should be regarded as regularised 
quantities. For a low-pass filter whose kernel decreases properly to zero for 
large values of its argument, Leray established uniqueness and regularity of 
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the solution u of the Leray-regularised equations, as well as convergence to a 
Navier-Stokes solution as the filter width tended to zero. This regularisation 
principle maintains the conservative structure of the filtered equations and 
represents a mathematically-based approach for obtaining reduced descrip- 
tions of turbulent flow. As such, Leray's regularisation principle may bridge 
the gap between practical requirements posed by LES and the high level of 
mathematical rigour required to guarantee systematic progress in turbulence 
analysis. The main practical question, of course, is whether or not this coars- 
ening of the flow description leads to accurate predictions of the smoothed 
flow. One of the goals of the present paper is to resolve this issue favourably 
in the context of turbulent mixing in a shear layer. 

In comparison to Leray's regularisation principle, the framework of the La- 
grangian averaged Navier-Stokes— a (LANS— a) equations provides an alter- 
native systematic method for modelling the mean circulatory effects of small- 
scale turbulence, while maintaining the mathematical properties that guar- 
antee existence of a unique, regular solution and a finite dimensional global 
attractor |13[ I14|, The inviscid Lagrangian averaged Euler— a equations were 
originally derived as Euler-Poincare equations in the framework of Hamil- 
ton's principle for geometric fluid mechanics JH]- The corresponding LANS— a 
model for turbulent flow may also be obtained more directly, by applying 
Lagrangian filtering to the Kelvin circulation theorem for the incompress- 
ible motion of a fluid with viscosity. As a consequence, the equations for the 
LANS— a model are obtained, and these generalise the Leray regularisation. 
The LANS— a solutions also converge to the unfiltered Navier-Stokes solutions 
as the filter width tends to zero. However, because the filtered Kelvin theorem 
is built into this turbulence modelling approach, it is also suitable for mod- 
elling turbulent flows in a rotating frame of reference and with buoyancy, as 
required in weather and climate modelling. The present paper compares the 
performance of the LANS— a model with the Leray model in numerical simula- 
tions of turbulent mixing and investigates the improvements in its predictions 
relative to the required computational effort. These findings are compared to 
traditional subgrid models in simulations obtained using the popular dynamic 
eddy- viscosity model [T5| . 

1.3 Outline of the paper 

This paper is organised as follows. Section [2] presents and analyses regulari- 
sation modelling for large eddy simulation of turbulent flow. The Leray and 
LANS— a models will be derived and discussed. Section |3] describes the tur- 
bulent mixing layer and the numerical methods adopted. The main features 
of the developing transitional and turbulent flow as captured by the regular- 
isation and dynamic subgrid models are presented to characterise the global 
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trends. Section 0] focuses on the quantitative predictions of turbulent mixing 
obtained with the Leray model. Mean flow, fluctuating quantities and spectral 
properties of the flow will be considered in order to assess the accuracy of 
the Leray model in some detail, both at low and at high Reynolds numbers. 
Section [5] is dedicated to improvements in the predictions that arise from the 
a-extensions of the Leray principle. It is shown that small-scale variability is 
restored to some degree in the LANS— a model, which goes at the expense 
of more strict requirements on the spatial resolution. Concluding remarks are 
collected in section 



2 Regular isat ion modelling for large eddy simulation 

In this section we relate mathematical regularisation of the convective fluxes 
to implied subgrid models for LES. First, in subsection 12.11 we review the 
filtering approach to large eddy simulation. Moreover, we describe the com- 
mon phenomenological treatment of the closure problem in which dissipative 
- or dispersive features of the turbulent stresses are represented by (dynamic) 
eddy-viscosity or similarity modelling respectively. This provides the context 
for discussing the basic mathematical regularisation strategy. The Leray and 
LANS— a regularisation principles provide central examples of this approach 
and the corresponding 'implied' subgrid models are discussed in subsection 
12.31 Finally, a simple Fourier-mode analysis of the regularisation models is 
collected in subsection 12.41 and compared with results obtained for standard 
similarity models. 

2.1 Filtering approach to large eddy simulation 

Filtering the Navier-Stokes equations requires a low-pass spatial filter L. Often, 
a convolution filter is adopted which in one spatial dimension associates the 
filtered velocity u with the unfiltered velocity u through 



with normalised filter-kernel G(z,£). This filter-kernel is characterised by an 
externally specified length-scale parameter i which defines the effective filter- 
width A, e.g., as PH 




(1) 
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where H(k, £) is the Fourier-transform of the filter-kernel. This definition ap- 
plies to all kernels that are square integrable such as the popular top-hat filter, 
the spectral cut-off filter or the Gaussian filter ^H]. Other definitions proposed 
in literature (see |Hj for an overview) are more restricted in their applicability 
to different filters. 

For incompressible fluids, the application of the filter L to the continuity 
and Navier-Stokes equations leads to: 

djUj = , and dtUi + dj(ujlli) + dip — T^^jj^i = —Oj(lIiUj — UjUi) (3) 

Here dt (resp. dj) denotes partial differentiation with respect to time t (resp. 
spatial coordinate Xj). Summation over repeated indices is implied. The com- 
ponent of the filtered velocity in the Xj direction is Hj and p is the filtered 
pressure. Finally, Re denotes the Reynolds number of the flow. 

In formulation (|3J) of the LES-template we recognise the application of 
the Navier-Stokes operator to the filtered solution {uj,p} on the left-hand 
side. On the right-hand side, the terms expressing the central closure problem 
are collected in the divergence of the turbulent stress tensor 
This tensor cannot be calculated from the filtered solution alone. Hence, one of 
the aims in the development of large eddy simulation is the effective capturing 
of the primary dynamical effects of Tij in terms of model tensors that may be 
evaluated in terms of operations on the filtered solution alone. 

In the absence of a comprehensive theory of how the small scales of tur- 
bulence influence its large scales, empirical knowledge about modelling Tij is 
essential, but it is still rather incomplete. Usually, LES subgrid models are pro- 
posed either on the basis of their presumed dissipative nature, or in view of the 
scale-similarity property of Tij in an inertial range ( 19 ). As further guidance 
in the construction of suitable models, one may attempt to incorporate in- 
formation associated with mathematical properties of the modelling problem 
such as realisability conditions ( 4 ), algebraic identities ( 5 ) or approximate 
inversion of the filter ( j^j, 20 ]). While realisability conditions may impose 
bounds on certain model parameters, the incorporation of algebraic proper- 
ties such as Germano's identity, possibly combined with filter-inversion 0, 
has led to a successful class of dynamic subgrid models. These dynamic 
models represent the current state of the art, at least for flows away from solid 
boundaries. 

Much of turbulence phenomenology is captured in the Kolmogorov picture, 
in which kinetic energy cascades in an average sense through an inertial 
range toward ever smaller scales, until its flow features are sufficiently lo- 
calised that they may be effectively dissipated by viscosity |21I22| , The dynam- 
ics is dominated by viscosity for flow features whose sizes (length scales) are 
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smaller than the Kolmogorov length t}k- When a filter of filter- width A S> r]K 
is applied to a turbulent solution, (virtually) all molecular dissipation associ- 
ated with the small scales is removed. In LES, the filter width is commonly 
chosen to be within a presumed inertial range in which the dynamics is domi- 
nated by convection. The continuous cascading of energy through the inertial 
range is then usually balanced by an "extra" eddy- viscosity contribution. This 
may yield a computational model that at least retains accurate large-scale in- 
formation in the filtered solution. Dissipation of turbulent kinetic energy was 
first parameterised by using the Smagorinsky model 9 in which one ap- 
proximates: 

nj - mg = -(C s A) 2 \S\S tJ (4) 

where Cs denotes Smagorinsky 's constant, Sij = dillj + djUi is the rate of 
strain tensor and IS*) 2 = SijSij is its magnitude. 

Besides the modelling of dissipation, a large volume of literature on the 
phenomenological treatment of the turbulent stress tensor is based on the 
assumed scale- similarity properties of this tensor when the filter-width 
is taken to be somewhere within the inertial range |19j . More precisely, this 
suggests that approximate models for can be obtained by applying its 
definition to an appropriate operation on the filtered solution. We may express 
as the commutator of the spatial filter and the product operator [HIE!, i.e., 

Tij = ujuj - UiUj = L(n„(u)) - Uij(L(u)) = [L, IIi 3 -](u) (5) 

where the product operator Ily(u) = u,- L Uj. The similarity aspects of the clo- 
sure problem in the inertial range were first parameterised by the Bardina 
model 23] ■ In this model one puts 

Tij -^mfj= UiUj - UiUj = [L, IIy](u) (6) 

Thus, the Bardina model applies the definition of Tij to the available filtered 
velocity. 

Extensions have been proposed in which the filtered velocity u is replaced by 
an approximate reconstruction of the unfiltered solution. In fact, for so-called 
graded filters L (e.g., including the top-hat and Gaussian filters, but excluding 
the spectral cut-off filter), operating on u with a formal approximation of the 
inverse L , leads to a partial reconstruction of the unfiltered solution, at least 
where (most of) the resolved scales are concerned. This reconstructed solution 
may be used to define generalised similarity models 0. This approach forms 
the basis for the approximate de-convolution models developed in |2Uj . These 
(generalised) similarity models characterise the dispersive effects of the small- 
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scale turbulence on the resolved flow features. When smooth test-solutions are 
considered and high-order approximate inversion, it may be shown that the 
difference between and the generalised model tensor scales as A m [Slll8j 
where the power m is controlled by the quality of the inversion. Of course, 
this scaling only holds in the mathematical limit A — ► 0; actual simulations 
with these models indicate the need for additional smoothing, e.g., through 
extra eddy-viscosity contributions. This is usually achieved in the context of 
the dynamic procedure to which we turn next. 

2.2 Dynamic mixed models 

The Smagorinsky and similarity models separately describe important intu- 
itive features of the turbulent stresses. However, these models are well known 
to be seriously flawed in their own ways. The Smagorinsky model displays 
low levels of correlation with and often leads to excessive dissipation, es- 
pecially near solid walls and in laminar flows with large gradients. This may 
even hinder a modelled flow from going through a complete transition to tur- 
bulence jllj . The similarity model of Bardina is known to display high cor- 
relation ^5], but it fails to provide effective dissipation of energy and it may 
give rise to unrealistically high levels of small scale fluctuations in the solu- 
tion. For these reasons, so-called mixed models have been proposed which 
combine similarity with eddy- viscosity models |24| . As an example, a basic 
mixed model is Tjj — > mM = mfj — CdA 2 \S\Sij in which denotes the eddy 
coefficient. The central problem that now arises is how this eddy coefficient 
should be specified in accordance with the evolving flow. A well-known and 
elegant way to approach this without unduly introducing ad hoc parameters 
is based on Germano's identity |16j which provides the basis for the dynamic 
subgrid modelling procedure. 

In dynamic models, the eddy-viscosity is intended to reflect local instan- 
taneous turbulence levels. One starts from Germano's identity: 

Tij — Tij = Rij , (7) 

where 

= U{Uj — TtiUj and Rij = (uiuj) — u % Uj . (8) 

Here, in addition to the basic LES-filter (•) of width A, a so-called explicit test 
filter (•) of width A is introduced. The only external parameter to be specified 
is the ratio of filter widths, which is often assigned as k = A/ A = 2. The 
implementation of the dynamic procedure starts by assuming a (mixed) model 



10 



Leray and LANS— a modelling of turbulent mixing 



rriij = dij(u) + cbij(u) for and Mj,- = Aij + CBij for Tij where Ay = a%j (u), 
5jj = 6jj(u). Here a^- and bij express assumed basic models, e.g., similarity 
or eddy-viscosity models, and the assumption C ~ c is usually invoked. If we 
introduce Aij = A^ — a^j, Bij = Bij — bij and use the common approximation 

(cbij) c(bij), insertion in Germano's identity yields Aij + cBij = Rij. This 
relation should hold for all tensor-components, which of course is not possible 
for a single scalar coefficient field c. To resolve this situation we introduce a 
linear averaging operator (•) and define the Germano residual by 



The averaging can involve integration over parts of the flow domain , over 
past history of the solution or include some Lagrangian averaging We 
hence obtain an optimality condition for c by requiring the variation of e to 
vanish. We can solve the local coefficient as |2(i| 



where we assumed (cfg) ~ c(fg). In order to prevent numerical instability 
caused by negative values of the eddy- viscosity, the dynamic coefficient is also 
artificially set to zero wherever (I10H returns negative values. This is referred 
to as 'clipping'. For further details we refer to |271IT%] . 

Dynamic subgrid models have contributed in many ways to the understand- 
ing of turbulent flow in complex situations. However, dynamic models are also 
known to be hampered by a number of drawbacks. First, the dynamic pro- 
cedure is quite expensive in view of the relatively large number of additional 
explicit filter operations that need to be included. Moreover, the implementa- 
tion contains various ad hoc elements or inaccurate assumptions such as the 
independence of the dynamic coefficient on the filter-level, or the well-known 
'clipping' of negative eddy-viscosity, required to ensure stability of a simula- 
tion. The achieved accuracy in actual simulations remains quite limited, e.g., 
due to shortcomings in the assumed base models. Since the dynamic approach 
does not contain external parameters other than the ratio between the width 
of the explicit test-filter and the basic LES filter, there is no chance of improv- 
ing the predictions by 'tinkering' with parameters. Finally, an extension to 
flows involving complex physics and/or developing in complex flow-domains 
is difficult since no systematic framework exists for this purpose. For these 
reasons, an alternative modelling approach is summoned and we turn to the 
recently proposed regularisation modelling ( j^Hj) in the next subsection. 





c = 




(10) 



Leray and LANS— a modelling of turbulent mixing 



11 



2.3 Regularisation strategy to subgrid closure 

In this subsection we consider two regularisation principles for the Navier- 
Stokes equations and derive the associated subgrid models in case the basic 
filter L has a formal inverse L~ l |28| . We consider Leray regularisation 
|T2] and the Lagrangian averaged Navier-Stokes-a (LANS— a ) approach 

The mathematical regularisation of the Navier-Stokes equations which we 
pursue here involves a direct and explicit alteration of the nonlinear connec- 
tive terms. In the context of this paper, this provides a systematic frame- 
work for deriving a subgrid model which is in sharp contrast with traditional 
phenomenological subgrid modelling. As sketched in the previous subsection, 
such phenomenological subgrid modelling achieves the desired smoothing of a 
turbulent flow only indirectly by quite 'independently' introducing a subgrid 
model for the turbulent stresses. 

Phenomenological subgrid models are usually only loosely connected to the 
basic LES filter L that was used to arrive at the unclosed equations. For ex- 
ample, eddy-viscosity models at best explicitly incorporate the width of the 
filter, and the modelled system does not contain any further information that 
characterises the shape of the particular filter. For similarity models there is 
a more direct relation with the explicit filter, although properties of the basic 
LES filter are never actually introduced, other than that the filter needs to be 
of convolution type. This 'independence' of the specific filter also holds for the 
dynamic procedure which only requires explicit specification of the test-filter. 
So, although the filter L is the only element that defines the relation be- 
tween DNS and LES, as well as the detailed properties of the subgrid stresses, 
the actual phenomenological modelling does not reflect this central role of 
L. Moreover, introducing dissipation by eddy-viscosity does not do justice to 
the intricacies of turbulent transport phenomena; it can at best perhaps char- 
acterise effective statistical properties of the kinetic energy dynamics. These 
fundamental considerations are fully respected in the regularisation principles 
that are considered next. 

Approach. A mathematical modelling approach for large eddy simulation can 
be obtained by combining a regularisation principle with an explicit filter 
and its formal inversion |28j . Historically, the first example of a smoothed 
flow description in this category is the Leray regularisation . Although this 
regularisation was introduced for entirely different reasons, we may reinterpret 
the Leray proposal in terms of its implied subgrid-model. This provides a direct 
connection with large eddy simulation. 

The application of a specific filter in combination with a mathematical reg- 
ularisation of Navier-Stokes equations allows one to incorporate properties 
from first principles into the modelled equations. Appropriate mathematical 
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regularisation results in a system of equations that is guaranteed to have a 
unique solution with suitable smoothness. The dynamics are characterised by 
an attractor of finite dimensions which has very favourable consequences for 
the computability of the solution, compared to the unfiltered Navier-Stokes 
system. Moreover, the momentum-conservation structure of the equations is 
retained. 



2.3.1 Leray regularisation. In Leray regularisation, one alters the connec- 
tive fluxes into TtjdjUi, i.e., the solution u is convected with a smoothed veloc- 
ity u. Consequently, the nonlinear effects are reduced by an amount governed 
by the smoothing properties of the filter operation, L. The governing Leray 
equations are [T^lHl] 

djUj = ; d t Ui + UjdjUi + dip - — Aui = (11) 

He 

Leray solutions possess global existence and uniqueness with proper smooth- 
ness and boundedness, whose demonstration depends on the balance equa- 
tion for J|u| 2 d s x. Based on the Leray equations (jll|) we may eliminate u by 
assuming u = L(u) and u = L _1 (u). For convolution filters one has, e.g., 
dtUi = d^L^iui)) = L^ 1 (dtUi) and the nonlinear terms can be written as 
Ujdj(ui) = djiiijUi) = dj (Hj L _1 . Consequently, one may readily obtain: 

L _1 (d t Ui + dj(ujlli) + dip - r^Aui^J = -dj (ujL~ l (tli) - L _1 (%Ui)^ (12) 

This may be recast in terms of the LES template as: 

d t Ui + dj(ujUi) + dip - jfe&Ui = -dj (mf^ (13) 

The implied asymmetric Leray model mfj involves both L and its inverse and 
may be expressed as: 

m\j = L (llj L~ x (ui) \ — Ujlli = UjUi — Ujlli (14) 

The reconstructed solution Uj is found from any formal or approximate inver- 
sion L~ l . For this purpose one may use a number of methods, e.g., polynomial 
inversion geometric series expansions |2Uj or exact numerical inversion of 
Simpson top-hat filtering [7]. 
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2.3.2 LANS— a regularisation by Kelvin filtering. A regularisation princi- 
ple which additionally possesses correct circulation properties may be obtained 
by starting from the following Kelvin theorem: 



where T(u) is a closed fluid loop moving with the Eulerian velocity u. The 
unfiltered Navier-Stokes equations may be derived from (|15[) |13U18j . This pro- 
vides some of the inspiration to arrive at an alternative regularisation principle 
for Navier-Stokes turbulence |13j . In fact, the basic regularisation principle was 
originally derived by applying Taylor's hypothesis of frozen-in turbulence 
in a Lagrangian averaging framework 29 . In this framework, the fluid loop is 
considered to move with the smoothed transport velocity u, although the cir- 
culation velocity is still the unsmoothed velocity, u. That is, in (|15j) we replace 
r(u) by r(u); so the material loop T moves with the filtered transport veloc- 
ity. From this filtered Kelvin principle, we may obtain the Euler-Poincare 
equations governing the smoothed solenoidal fluid dynamics, with djllj = 
and 



Comparison with the Leray regularisation principle in (|11|) reveals two addi- 
tional terms in ()16l) . These terms guarantee the regularised flow to be consis- 
tent with the modified Kelvin circulation theorem in which T(u) — > r(u). For 
LANS— a the analytical properties of the regularised solution are based on the 
energy balance for f u • L(u) d 3 x. 

The Euler-Poincare equations (|16|) can also be rewritten in the form of the 
LES template. The extra terms that arise in (|16[) give rise to additional terms 
in the implied subgrid model: 



d t Ui+dj(ujUi)+dip-—Aui = -dj{ujUi-UjUi)--(ujdiUj - UjdiuA (17) 



We observe that the Leray model (|14j) reappears as part of the implied 
LANS— a subgrid model on the right-hand side of (|17jl . Compared to the 
Leray model, the additional second term in the LANS— a model takes care 
of recovering the Kelvin circulation theorem for the smoothed solution. This 
formulation is given in terms of a general filter L and its inverse. After some 
further rewriting it may be shown that this model can be formulated in con- 
servative form, i.e., a tensor mfj can be found such that the right hand side of 
(fT7|) can be written as —djmf;. We illustrate this next for a particular filter. 




(15) 



d t Uj + u k d k Uj + UkdjUk + djp - dj(-u k u k ) - Tr^uj = (16) 



1 
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The subgrid model presented in (|17|) can be specified further in case the filter 
L has the Helmholtz operator as its inverse, i.e., Ui = L~ l (lli) = (l—a 2 djj)Hi = 
He a (uj). Then we recover the original LANS— a equations [Ej- The LANS— a 
model derives its name from the length-scale parameter a ~ A/5 [HOI- After 
some rewriting, the following parameterisation for the turbulent stress tensor 
is obtained |3"T] : 



The first term on the right-hand side is the Helmholtz- filtered tensor-diffusivity 
model ^Sj. The second term combined with the first term, corresponds to 
Leray regularisation using Helmholtz inversion as filter. The third term com- 
pletes the LANS— a model and maintains Kelvin's circulation theorem. In (|18jl 
an inversion of the Helmholtz operator He a is required which implies appli- 
cation of the exponential filter |32j . However, since the Taylor expansion of 
the exponential filter is identical at quadratic order to that of the top-hat or 
the Gaussian filters, one may approximate He" 1 , e.g., by an application of the 
explicit top-hat filter, for reasons of computational efficiency |3Uj . 

Relation to LES. Although the regularised turbulence equations (|13|) and 
(|17|) are formally similar to LES equations, they arose from different prin- 
ciples, as sketched above. Through the combination of an explicit filter and 
its inversion, the regularisation principle allows a systematic derivation of the 
implied subgrid-model. This resolves the closure problem consistent with the 
adopted filter. Even though the Leray and LANS— a formulations may, in ret- 
rospect, be interpreted in terms of implied subgrid- models, the regularisation 
equations were not obtained by applying the LES filtering method. Instead, 
this approach to modelling turbulence from the viewpoint of mathematical 
regularisation is an alternative to LES. 

In the next subsection we will adopt simple Fourier-mode analysis in one 
spatial dimension, in order to illustrate some basic properties of the regulari- 
sation models, before assessing these models in actual large eddy simulations 
of turbulent mixing in the next sections. 



2.4 Fourier- analysis of regularisation models 

In order to illustrate the dynamic effects of some of the subgrid models intro- 
duced in either of the previous subsections, we will first investigate the subgrid 
flux due to a single Fourier-mode. Specifically, we select u = sin(A;x) which im- 
plies u = H(kA)u for convolution filters whose kernel G has Fourier-transform 




(18) 
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H. Substitution in the definition of the turbulent stress yields 



r = u 2 — w 




The corresponding flux f T 



d x r may be expressed as: 



f T = k[H 2 (kA) - H(2kA)) sm(2kx) = kA T (kA) sin(2kx) (20) 



in which we introduced the characteristic amplitude function A T which de- 
pends on kA only. We return to this amplitude function momentarily and use 
it to assess different subgrid models. 

The various subgrid models introduced above will now be evaluated in the 
same manner as in (|19|) and (|2Uj) . In one spatial dimension the Leray and 
LANS— a models yield an identical flux when applied to u = sin(kx). For the 
Leray model one may show after some manipulation 



f L = k(H 2 (kA) - H(kA)H{2kA)) sin{2kx) = kA L (kA) sin(2/cx) (22) 



Comparing (|22j) with ()20|) we observe a close resemblance with the only differ- 
ence due to a term H(kA) arising where unity appears in the exact expression. 
Since the filter is assumed normalised, we have H(0) = 1 and hence the devi- 
ations will be small provided A; A <C 1. For larger values of kA the deviations 
may be more significant. We return to these points shortly and first complete 
the expressions for the other subgrid models. 

For a single Fourier mode the Bardina similarity model results in 





with corresponding flux 



tub = u 2 — u = 




with corresponding flux 



f B = kH 2 (kA)(H 2 (kA) - H(2kA) \ sin(2fcx) = kA B (kA) sm(2kx) (24) 
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We observe that in this single Fourier-mode evaluation mg = H 2 (kA)r, show- 
ing close agreement as kA <C 1. In three spatial dimensions, the Bardina model 
requires a number of additional applications of the explicit filter, which adds 
considerably to its computational cost. An approximation to the Bardina sim- 
ilarity model is the tensor-diffusivity model [33] which does not require the 
extra explicit filtering. We consider this model in more detail next. 

The tensor-diffusivity model may be derived by truncating a formal Tay- 
lor expansion of the turbulent stress tensor |34| . We express a general one- 
dimensional integral filter as: 

/(*)=/ )/(£)#= / G(n)f(x + A V )d V (25) 

•J — oo J — OO 

where the second expression is in terms of r/ defined by Arj = £ — x. We 
assume that the signal / has a globally convergent Taylor expansion so that 
for normalised filter- kernels the application of the filter may be expressed as: 

/ = / + V A m M m f^ ; M m = — G( V ) V m dr, (26) 

where /( m ) denotes the m-th derivative of / and we introduced the moments 
M m of the filter, including a factor 1/ml for notational convenience. This 
expresses the filtered signal in terms of derivatives of the unfiltered signal. 
Application to the turbulent stress tensor in one spatial dimension allows to 
write 

r = u 2 — u 2 

oo oo 2 

= {u 2 + Y, A m M m (u 2 )( m >} ~{u+J2 & m M m u {m) } 

771=1 771=1 

OO OO OO 

= ^ A m M m {(u 2 ) {m) - 2uu^} - Y,Y, Am+nMmMnu(m)uin) ( 27 ) 

771=2 771=1 77=1 

where use was made of (u 2 )' — 2uu' = 0. This expression for the turbulent 
stress holds for general filters for which all moments M m exist. As an example, 
this includes popular filters such as the top-hat, Helmholtz or Gaussian filter 
[TS] . but it is equally valid for general compact support high-order filters [HJ. 
Collecting terms according to their power of A we have after some rewriting 

r = A 2 (ti') 2 {2M 2 - M x 2 } + A 3 (V) 2 )'{3M3 - MiM 2 } 
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+ A 4 j (V) 2 )' 4M 4 - M1M3 +(u"f 2MiM 3 - M| - 2M 4 } + . . .(28) 

This expression contains contributions from all terms ~ A n with n > 2 and 
involves both even and odd order moments. Many popular filters are formu- 
lated in terms of an even filter- kernel, i.e., G{—rj) = G(rj). For these filters the 
odd moments are identically zero and correspondingly all contributions ~ A" 
with n odd, are zero as well. In case the filter is assumed to be an iV-th order 
filter, i.e., by definition M m = 5 m o for m = 0, 1, . . . ,N — 1 at some N > 1, a 
corresponding additional number of contributions of lower order in A is equal 
to zero. This provides some control over the formal magnitude of the turbu- 
lent stress tensor and may be interpreted as a less strong filtering in the sense 
that the effective filter-width defined in ((2j) decreases rapidly with increasing 
filter-order |17j . 

For first - or second order filters that are commonly considered in large eddy 
simulation, the truncation of (|28|) at order A 2 yields the standard tensor- 
diffusivity model. To replace derivatives of u in this expression with corre- 
sponding derivatives of u it is required to approximately invert the filter. Con- 
sistent to the required order in A, we may approximate the inversion of the 
filter simply as u ~ u which leads to 

m T D = A 2 {2M 2 - Mf}(d x ll) 2 = A 2 V 2 (d x ll) 2 

= ^{2M 2 - M 2 }{kA) 2 (^H 2 (kA) + H 2 (kA)cos(2kx)^j (29) 

with corresponding flux 



f TD = k([2M 2 - Mf}(kA) 2 H 2 (kA)^j sin(2A;x) = kA TD {kA) sin(2A;x) (30) 

and, for notational convenience, we introduced the filter variance V 2 = 
2M 2 — M 2 . In case a higher order accurate consistent truncation of (|28|) is 
desired, or, when the filter itself is of higher order, then the expression of u in 
terms of u and its derivatives becomes slightly more involved. However, this 
poses no principal problems and one may invert the class of filters consid- 
ered here at least formally up to arbitrary order. We will not consider these 
extensions but restrict ourselves to the standard tensor-diffusivity model in 
(Hi). 

The tensor-diffusivity model (|29|) is known to induce instabilities into a 
large eddy simulation [311131?] . These instabilities may be removed by adding a 
(dynamic) eddy-viscosity contribution. However, the instabilities may also be 
addressed by turning to the spatially filtered tensor-diffusivity model instead, 
as was observed in numerical experiments [30] and recently established analyt- 
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ically [35]. In fact, it was shown that spatial filtering of the tensor-diffusivity 
model provides a closure that yields uniqueness and global regularity. For the 
filtered tensor-diffusivity model we find: 



We notice the extra smoothing of the flux in (|32|) compared to (j30[) . arising 
from the extra filter which adds a factor H(2kA). For all popular large eddy 
filters this extra filter appears to damp the behaviour of the flux at high 
wavenumbers sufficiently to restore boundedness of the solution, however, at 
the expense of increasing the computational costs. 

In order to assess some of the properties of the regularisation and similarity 
subgrid models introduced above, we begin by considering the implications 
in case kA <C 1 and turn to properties at large kA momentarily. In the low 
wavenumber range the subgrid fluxes are small and we can obtain a precise 
impression of the type of deviations that arise. By Taylor expansion of the 
characteristic amplitude functions we have: 



m fTD = A 2 V 2 (d x u) 2 

= ^V 2 {kA) 2 ^H 2 (kA) + H 2 (kA)H{2kA) cos(2A;x)) (31) 



with corresponding flux 




A L (z) 



A T (z) 



( 



(H"(0) - H'(0) 2 y 2 - (H'"{0) - #'(0)F"(0))z 3 
L H ""{Q) - l -H'{Q)H"'{Q) - \H"(0) 2 )z 4 + ... 

■H'(0)z- (Jtf"(0)+#'(0) 2 )^ 2 



(33) 






(iT"(0) - ^'(O) 2 ^ 2 - (V'(0) + ff'O 
1 H ""(0) + ~H'(0)H">(0) + ^H"(0) 2 
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- 2H'(0) 2 H"(0) - #'(0) 4 ) z 4 + . . . (35) 
Atd(z) = V 2 z 2 + 2V 2 H'(0)z 3 + V 2 (V(0) 2 + H"(pj\ z 4 + . . . (36) 
A fTD (z) = V 2 z 2 + W 2 H'(0)z 3 + 1/ 2 (W(0) 2 + 3F"(0))z 4 + . . . (37) 

where we assumed that all filters are normalised, i.e., H(0) = 1. If we assume 
in addition that the filter-kernel G is an even function, as is common in most 
large eddy studies, then all odd order derivatives H ( 2n+1 ) (0) = and the above 



expansions simplify to 

M*) = ~H"(0)z 2 - (^ff""(0) - \h"(0) 2 )z 4 + . . . (38) 

A L {z) = - 3 -H"(0)z 2 - f^H""(0) + 3 -H"(0f)z 4 + ... (39) 

A B {z) = -H"(0)z 2 - (^2?""(0) + \H"{V) 2 )z 4 + ... (40) 

Atd(z) = V 2 z 2 + V 2 H"(0)z 4 + ... (41) 

A fTD (z) = V 2 z 2 + 3V 2 H"(0)z 4 + ... (42) 



The expansions of the turbulent stress tensor in (|33|) or ()38j) provide a point 
of reference for the different models. The Leray model gives rise to a term ~ z 
in (|34|) . which is absent in the full turbulent stress tensor. Also, the second or- 
der contribution in (|34|) or (|39|) deviates systematically from the corresponding 
term in (|33[) . Higher order corrections deviate as well but, since 7/12 ~ 5/8, 
corrections involving the fourth order derivative of H are still a fair approx- 
imation. These expansions establish that at low kA properties of the Leray 
(and LANS— a) model deviate markedly from properties of the exact closure 
problem. As mentioned before, this illustrates that the regularisation approach 
may be re-interpreted in the language of the spatial filtering formulation for 
LES but it may not be directly derived from it. 

The Bardina similarity model corresponds exactly at second order with 
the full turbulent stress tensor. Moreover, for a symmetric filter the fourth 
order corrections deviate from the exact turbulent stress in a way that is 
quite comparable to the Leray model. The expansions for the tensor-diffusivity 
model appear different at first sight. However, it is not hard to verify that 
V 2 = —(H"(0) — -fT'(O) 2 ) and hence we notice that the lowest order terms for 
both the tensor-diffusivity and the filtered tensor-diffusivity model are exactly 
the same as in the full turbulent stress. Higher order corrections deviate in a 
way similar to what is observed in the Bardina model. 
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The truncated expansions are meaningful only for sufficiently low wavenum- 
bers k. In this wavenumber range the subgrid flux is not particularly impor- 
tant for the dynamics of the flow. To further assess the suitability of a subgrid 
closure for LES, the properties in the high wavenumber range are more impor- 
tant. We consider these next by explicitly evaluating the dependence of the 
characteristic amplitude functions on kA for three well known symmetric LES 
filters, i.e., the top-hat (th), the Helmholtz (H) and the Gaussian (G) filters. 
Specifically, the Fourier transform of these filters is given by 

H th (z) = Sm ^ /2) ; H H (z) = \— ; H G (z) = exp(-— ) (43) 

The leading order expansion of each of these filters is given by 1 — z 2 /24 + . . . 
. However, these filters differ essentially in their attenuation of high-/c solution 
components. As z 3> 1 the amplitudes of these Fourier-transforms yield either 
a relatively slow algebraic decrease \Hth\(z) ~ z" 1 or Hjj(z) ~ z~ 2 , or a very 
fast exponential decrease for the Gaussian filter. 




Figure 1. Characteristic amplitudes as function of kA for the symmetric top-hat (a), Helmholtz 
(b) and Gaussian (c) niters: turbulent stress (solid), Leray (dashed), Bardina (dash-dotted), 
tensor-diffusivity (dotted) and filtered tensor-diffusivity (o). 



The properties of the filter directly translate into the behaviour of the char- 
acteristic amplitude functions. Correspondingly, also the subgrid flux associ- 
ated with the various models is influenced. To illustrate this, we plotted the 
amplitude functions versus kA for the three filters (j43j) in figure ^ As point 
of reference the characteristic amplitude of r is shown as solid line, indicating 
either an oscillatory behaviour associated with the top-hat filter or a uni- or 
bi-modal dependence on kA for the Gaussian and Helmholtz filters, respec- 
tively. The decay of the amplitude function for large kA differs in accordance 
with the asymptotic behaviour of the Fourier transforms; the Gaussian filter 
results in the smallest 'active' range of wavenumbers, while the top-hat filter 
provides the largest range. 
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The Leray model is seen to deviate already at quite low A; A, consistent with 
the Taylor expansion shown above. However, for the top-hat and Gaussian 
filters the correspondence of the Leray model with the exact turbulent stress 
is better and over a wider range than that observed for the other models. 
For the top-hat filter the maximum amplitude arises at slightly too low kA, 
indicating that the Leray model may exaggerate the larger scales in a solution. 
In combination with the Gaussian filter the tail of the amplitude function is 
nearly perfectly predicted by the Leray model, although there is a similar over- 
prediction of the maximum amplitude at slightly too low wavenumber. For the 
Helmholtz filter the correspondence of the Leray model is still quite close but 
the filtered tensor-diffusivity model agrees somewhat better with the exact 
amplitude function. None of the models properly captures the change of sign 
in the amplitude function of r that arises in combination with the Helmholtz 
filter. 

Turning to the Bardina model, we observe that the amplitude is strongly 
under-predicted and that the spectral support is too restricted, i.e., the values 
of kA for which the amplitude is sufficiently different from zero is too small. 
The tensor-diffusivity model overestimates the amplitude in combination with 
the Helmholtz or the Gaussian filters and shows a periodic dependence on 
kA in combination with the top-hat filter. In particular, this means that, in 
combination with the top-hat filter, the amplitude function does not reduce to 
zero even for very small length-scales. These properties correlate directly with 
the underlying unstable nature of this model in actual simulations |34U35| . 
The filtered tensor-diffusivity model improves essentially on this behaviour 
but, apart from the combination with the Helmholtz filter, never seems to 
achieve accurate levels of subgrid flux. 

Although this single-mode analysis cannot be fully conclusive regarding the 
behaviour of any of these subgrid models in actual large eddy simulations, it 
does provide some interesting first illustrations that will be seen to correlate 
well with the findings in turbulent mixing flow to which we turn in the next 
section. 



3 Dominant flow-features in turbulent mixing 

In this section we first describe the turbulent mixing flow that is used in this 
paper to assess the quality of regularisation modelling. In subsection 13.11 we 
sketch the numerical method and illustrate the global development of this 
canonical flow. Then, in subsection 13.21 we present the reference large eddy 
simulation results obtained with the Leray and LANS— a models, restricting 
ourselves to a discussion of the prominent dynamic flow-structures that arise 
and how well these are captured with these models. A quantitative analysis is 
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presented in sections 0] and 03 



3.1 Numerical simulation of temporal mixing 

The governing equations are discretise using the so-called method of lines. 
In the computations we consider the compressible formulation and perform 
simulations at a low convective Mach number which was shown to provide 
essentially incompressible flow-dynamics We write the Navier-Stokes or 
LES equations concisely as dtU = TiJA) where U denotes the state- vector con- 
taining, e.g., (filtered) velocity and pressure, and T is the total flux, composed 
of the convective, the viscous, and in case of LES also the subgrid fluxes. 

The operator T contains first and second order partial derivatives with 
respect to the spatial coordinates. The equations are discretised on a uni- 
form rectangular grid and the grid size in the Xj-direction is denoted by hj. 
If we adopt a spatial discretization around a grid point x^*., the operator 
J-(JA) can be approximated in a consistent manner by an algebraic expression 
Fijk({U a f3 y }) where {C/ Q/ g 7 } denotes the state vectors in all the grid-points 
that cover the domain. Usually, only neighbouring grid points around (i,j,k) 
appear explicitly in i 7 ^. After these steps we obtain a large system of ordinary 
differential equations 

= F ijk ({U a ^}) ; U ijk (0) = U® (44) 

where Uij^t) approximates £/(xy&, t) and uj^l represents the initial condition. 
Hence, in order to specify the numerical treatment, apart from the initial and 
boundary conditions, the spatial discretization which gives rise to Fijk and the 
temporal integration need to be defined. We next introduce these separately 
which is central to the method of lines. 

The time stepping method which we adopt is an explicit four-stage compact- 
storage Runge-Kutta method JH]- When we consider the scalar differential 
equation du/dt = f(u), this Runge-Kutta method performs within one time 
step of size St 



u 



W = „(0) +W( „O--l)) ( i = 1,2,3,4) (45) 



with u(°) = u(t) and u(t + St) = . With the coefficients ft = 1/4, ft = 1/3, 
ft = 1/2 and ft = 1 this yields a second-order accurate time integration 
method |36| . The time step is determined by the stability restriction of the 
numerical scheme. It depends on the grid-size h and the eigenvalues of the 
flux Jacobi matrix of the numerical flux /. In a short-hand notation one may 
write St = CFL h/\X max \ where |A maa; | denotes the eigenvalue of the flux 
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Jacobi matrix with maximal size, and CFL denotes the Courant-Friedrichs- 
Levy-number which depends on the specific choice of explicit time integration 
method. For the present four-stage Runge-Kutta method a maximum CFL 
number of 2.4 can be established using a Von Neumann stability analysis. In 
the actual simulations we use CFL = 1.5, which is suitable for both DNS and 
LES, irrespective of the specific subgrid model used. 

In order to specify the spatial discretization we distinguish between the 
treatment of the convective and the viscous fluxes. We will only specify the nu- 
merical approximation of the c?i-operator; the 82 and c^-operators are treated 
analogously. Subgrid-terms are discretised with the same method as the vis- 
cous terms. Throughout, we will use a second order finite volume method [SZ1 
for both the viscous and the convective fluxes. The discretization of the convec- 
tive terms is the cell vertex trapezoidal rule, which is a weighted second-order 
central difference. In vertex (i,j,k) the corresponding operator is denoted by 
d\ and for the approximation of d\f it is defined as 

(dif)i,j,k = (si+ij,k ~ Si-i )jik )/ (2hi) (46) 
with Sij )k = (9i,j-i,k + 2gi,j,k + 5 r i,i+i,fc)/4 
and gijfi = (fij t k-i + %fi,j,k + fi,j,k+l)/^ 

This illustrates the second-order central difference applied to a quantity Sijk 
which itself is a local average of / over the X2 and X3 directions. The viscous 
terms contain second-order derivatives which are treated by a consecutive 
application of two first order numerical derivatives. This requires for example 
that the gradient of the velocity is calculated in centres of grid-cells. In centre 
(* + 5>i + 3>^ + !) the corresponding discretization D\f has the form 

(^i/X'+fj+i.fc+i = (Si+ij + i.fc+| - a ij+±,k+±)/hi ( 47 ) 
with Sjj + i )fc+ i = (/ij,fc + fi,j+i,k + fi,j,k+i + /t,j+i,fc+i)/4 

The second derivative is subsequently calculated in the point (i,j,k) by ap- 
plying the operator D\ again, but now to the staggered approximations of 
the first derivative. Specifically, we determine which, according to 

(|47|). contains terms s i+ ij >fe . These may be evaluated from the now known 

approximations at (i + + h, k + 5) which are available after the first appli- 
cation of D\ as defined in l[47[) . This basic numerical method, i.e., consisting 
of explicit second order accurate four-stage Runge-Kutta time-stepping and 
second order finite volume discretization was adopted for the direct and large 
eddy simulations of a turbulent mixing flow to which we turn next. 
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Problem statement. The three-dimensional temporal mixing layer is consid- 
ered at a Reynolds number based on upper stream velocity and half the initial 
vorticity thickness of 50. This is sufficiently high to allow a mixing transition 
to small scales. It is also sufficiently low to enable an accurate DNS that re- 
solves all relevant turbulent scales on the computational mesh. The governing 
equations are solved in a cubic geometry of side £ which is set equal to four 
times the wavelength of the most unstable mode according to linear stabil- 
ity theory {t = 59 in the present units). Periodic boundary conditions are 
imposed in the streamwise {x\) and spanwise (x^) direction, while in the nor- 
mal (X2) direction the boundaries are free-slip walls. The initial condition is 
formed by mean profiles corresponding to constant pressure p, u\ = tanh(x2) 
for the streamwise velocity component and U2 = U3 = 0. Superimposed on 
the mean profile are two- and three-dimensional perturbation modes obtained 
from linear stability theory. Further details of the problem statement for the 
three-dimensional temporal mixing layer may be found in 

Visualisation of the DNS data demonstrates the roll-up of the fundamental 
linear instability and successive pairings. Four rollers with mainly negative 
spanwise vorticity emerge from the initial condition (t = 20). After the first 
pairing (t = 40) the flow has become highly three-dimensional. Another pairing 
(t = 80) yields a single roller in which the flow exhibits a complex structure 
with many regions of positive spanwise vorticity. This structure is a result of 
the transition to turbulence which has been triggered by the pairing process 
at t = 40. The simulation is stopped at t = 100, since the single roller at 
t = 80 does not undergo another pairing in this computational model. The 
accuracy of the simulation with 192 3 cells is satisfactory for our purposes as 
was demonstrated by comparing with computations using, e.g., a fourth order 
discretization, or coarser grids containing 64 3 and 128 3 cells |27| . 

In the next subsection we will further specify the computational modelling 
for the large eddy simulation of this turbulent mixing flow and present some 
results that characterise the global flow-features as predicted by the regulari- 
sation models. This provides a first, general assessment of the quality of these 
models. 

3.2 Capturing instantaneous mixing flow with regularised large eddy 
simulation 

Throughout this paper, we consider large eddy simulations on three different 
grids with 32 3 , 64 3 and 96 3 cells, while keeping the filter- width fixed at A = 
1/16. This value of A provides a challenging test-case for large eddy simulation, 
studied first in JTT]. Here, t denotes the length of the side of the cubical 
computational domain that was used. In this way the subgrid resolution 
r = A/h is varied explicitly, and we can separately influence the effects of 
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numerical errors in the description. Moreover, we may assess independently the 
quality with which the LES models capture the physical features of the flow by 
considering the predictions obtained for the (approximately) grid-independent 
situation that is obtained if r is large enough |38j . 

The possible numerical contamination of the smaller resolved scales was 
found to be well characterised by the subgrid resolution. In the simulations 
considered here r varies between 2 and 6; for the lower value of r the specific 
spatial discretization has some effect, even for mean flow properties, while 
for r ~ 4 — 6 a more acceptable solution is obtained that is approximately 
grid-independent for most intends and purposes. The deviations from filtered 
DNS data that remain at r = 6 are a direct measure of the deficiencies in the 
subgrid model only and would also be found in the same way for other spatial 
discretization methods. 

The regularisation models require explicit filtering and (approximate) inver- 
sion. Moreover, the dynamic procedure that will be included for comparison 
purposes, requires explicit (test-)filtering as well. Various, filters can be used. 
In this paper, we will adopt invertible numerical quadrature to approximate 
the top-hat filter. In this way we arrive at a consistent representation of the 
top-hat filter and at the same time guarantee exact numerical inversion on 
the computational grid. In one dimension, symmetric numerical convolution 
filtering u = G * u corresponds to kernels 

771 

G{z)= aj S(z-Zj) ; \zj\< A/2 (48) 

j=-m 

defined on 2m + 1 points. In particular, in this paper we consider the simplest 
three-point filters with ao = 1 — s, at = a_i = s/2 and zq = 0, Z\ = —Z-\ = 
A/2. Here we use s = 1/3 which corresponds to Simpson quadrature of the 
top-hat filter. In actual simulations the resolved fields are known only on a set 
of grid points {x m }^ =0 . The application of the inverse L _1 of this three-point 
filter can be specified using discrete Fourier transformation. In particular, for 
a general discrete filtered solution {u(x m )} we have [Z] 

j=-n 

where the subgrid resolution r = A/h is assumed to be even. An accurate 
and efficient inversion can be obtained with only a few terms, recovering the 
original signal to within machine accuracy with n ~ 10. Filtering and inversion 
in three dimensions arises from composing three one dimensional filters. 
A first qualitative test of subgrid models is obtained by studying the predic- 
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tion of the dominant flow-structures in instantaneous solutions. So far in liter- 
ature, failure of a large eddy simulation to predict these global features of an 
evolving flow was not considered a decisive stumbling block for the particular 
subgrid model that was adopted. Typically, it is argued that differences in an 
initial DNS-field concerning length-scales well below the filter-width A would 
not be of importance to the corresponding filtered initial condition for LES. 
However, sensitive dependence of a turbulent flow to small-scale differences in 
initial conditions could result in significant effects on the instantaneous solu- 
tions at later times. Hence, the correspondence between LES and DNS does 
not necessarily have to extend to instantaneous realizations of turbulent fields. 
This point is certainly valid. However, at the modest Reynolds number consid- 
ered here, and with the laminar initial condition that is used, one may expect 
that accurate subgrid models will yield a strong correlation between filtered 
snapshots of the DNS solution and large eddy predictions. Correspondingly, a 
relevant and actually quite severe test-case for LES, albeit at modest Reynolds 
numbers is obtained. 

As a typical illustration of the mixing layer, the filtered DNS prediction of 
the vertical velocity and the corresponding Leray and LANS— a results can be 
compared in snapshots. We approximately eliminated the spatial discretization 
effects by using a resolution of 96 3 in both large eddy simulations. We observe 
that even at the instantaneous solution level, both the Leray and LANS— a 
models capture the 'character' of the filtered solution. While the Leray solution 
appears to provide a slight under-prediction of the influence of the small scales, 
the LANS— a solution corresponds more closely with the filtered DNS findings, 
restoring some of the small scale variability. These instantaneous predictions 
are both much better than those obtained with the dynamic eddy-viscosity 
model which turn out to be too smooth. 

The different regularisation models are known to have different effects on the 
tail of the resolved kinetic energy spectrum E(k). In the Kolmogorov picture 
of homogeneous, isotropic turbulence an inertial range in which E(k) ~ A;~ 5 / 3 
develops over an extended range of wavenumbers k up to a Kolmogorov 
wavenumber k v ~ I/77 where 77 is the viscous dissipation length-scale. This 
entire dynamic range needs to be properly captured in order to arrive at a re- 
liable DNS. The Leray and LANS— a models give rise to a spectrum in which 
there is a smooth transition from a —5/3 power law to a much steeper al- 
gebraic decay, beyond wavenumbers ~ 1/A. The sharper decrease of kinetic 
energy with wavenumber implies a corresponding strong reduction in required 
computational effort needed for the simulation of the relevant dynamic range. 
The LANS— a model displays a tail of the spectrum ~ k~ 3 while the Leray 
model decays even more steeply, as ~ 

fc-13/3 [T3)rra|. The steeper decay using 
the Leray model is directly reflected in the smoother impression of instanta- 
neous solutions. Hence, through the selection of A a direct external control is 
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achieved over the computational costs associated with the regularisation mod- 
els. This is illustrated in figure |5J In case an energy range of, say, m decades 
is desired then all wavenumbers up to fcz,(m), k a (m) and kDNs{fn) need to be 
resolved for the Leray, LANS— a and DNS approaches respectively. This cor- 
responds to a significant difference in the associated computational expense, 
while all three simulations would provide excellent accuracy at least for all 
wavenumbers up to ~ 1/A. 
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Figure 2. Sketch of resolved kinetic energy spectrum in a homogeneous, isotropic turbulence, 
displaying a —5/3 tail in DNS (solid), a —3 tail in LES using the LANS— a model (dashed) and a 
— 13/3 tail in LES using the Leray model (dash-dotted). 

In the next section we will turn to a more quantitative assessment of the 
Leray model and investigate a range of flow features including mean flow as 
well as spectral properties. 

4 Leray predictions of turbulent mixing 

In this section the Leray predictions for the turbulent mixing layer are dis- 
cussed. First, in subsection 14.11 we restrict ourselves to a modest Reynolds 
number and compare results directly with the available filtered DNS data. 
Subsequently, in subsection 14.21 we turn to much higher Reynolds numbers 
and establish the robustness and reliability of this computational model. 

4.1 Mixing at modest Reynolds numbers 

In order to assess the quality of a subgrid model it is important to address 
the accuracy of large eddy predictions for a variety of flow properties that 
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characterise different length-scale ranges. Moreover, a clear distinction should 
be made between the quality of the subgrid model and effects arising from the 
nonlinear interaction with discretization errors at marginal spatial resolution 
|38ll39j . We comply with these two general requirements by including LES 
predictions ranging from mean flow quantities and statistics of fluctuations to 
spectra of turbulent kinetic energy, at different spatial resolutions. Moreover, 
we will consider the 'inner operations' of the Leray model by studying the 
contributions of the subgrid flux to the energy dynamics and establish the 
amount of forward and backward scatter. 

4.1.1 Mean flow properties.. The resolved kinetic energy E and the mo- 
mentum thickness of the resolved flow 8 are two characteristic mean flow prop- 
erties of the temporal mixing layer. These quantities are defined as 

f 1 1 f^fo 

E= -UiUid* ■ 8=- / (l-(ui)(x2,*))((«i)(a?2,t) + l)d!C2 (50) 

where fi is the flow domain and (•) denotes averaging over the homogeneous 
xi and X3 directions. The evolution of E illustrates the transitional flow and 
subsequent self-similar decay in the turbulent regime. It depends primarily on 
the larger scales in the flow. Similarly, the momentum thickness is a large-scale 
quantity that can be used as a measure for the progress of the mixing. 




Figure 3. Resolved kinetic energy E (a) and momentum thickness <5 (b): filtered DNS (o), 
Leray- model (32 3 : dash-dotted, 64 3 : solid, 96 3 : A), dynamic model (32 3 : dashed, 64 3 : dashed with 
o). A fixed filterwidth of 1/1Q was used at a Reynolds number Re = 50. 



The simulation results for E and 5 are collected in figure 01 In an earlier 
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LES study of this flow by Vreman et al. it was established that the dy- 
namic eddy-viscosity model was among the subgrid models that performed 
best, compared to similarity - and other dynamic (mixed) models. We observe 
that the kinetic energy E is quite well represented during the laminar stages, 
but deviates considerably from the filtered DNS data during the turbulent 
stages. In particular, we notice that this deviation becomes larger in case the 
subgrid resolution is increased at fixed filter- width. This indicates that the 
grid-independent solution corresponding to the dynamic eddy-viscosity model 
does not provide sufficient dissipation on its own. Instead, the interaction 
with the particular spatial discretization method on coarser grids is seen to 
provide additional decay of E. The capturing of the momentum thickness with 
the dynamic model is not very accurate in the important self-similar turbu- 
lent regime: the grid-independent dynamic model yields a significant under- 
prediction of the mixing rate. 

The Leray results indicate that the resolved kinetic energy is also overesti- 
mated in the turbulent regime. For the grid-independent solution the predicted 
decay rate dE/dt is nearly constant in the turbulent stages and corresponds 
well with the filtered DNS data. Moreover, the interaction with the spatial 
discretization method on coarse grids is seen to result in an increased decay 
rate, similar to what was observed for the dynamic model. The Leray predic- 
tion for the momentum thickness compares significantly better with filtered 
DNS results than was obtained with the dynamic model. We observe that 
some of the discrepancies between Leray and filtered DNS results are due to 
numerical contamination on coarse grids. By increasing the resolution at fixed 
filterwidth A, a good impression of the grid-independent solution to the mod- 
elled equations can be inferred using 64 3 - 96 3 grid-cells, i.e., A/h = 4 to 
6 |H8j . Numerical contamination also plays a role in the dynamic model. The 
grid-independent prediction for 5 corresponding to the dynamic model appears 
less accurate than the corresponding Leray result. 

4.1.2 Statistics of velocity variations.. In order to assess the quality of 
the Leray predictions for quantities that are more sensitive to small-scale in- 
formation we next consider velocity variations defined by: 



where (•) denotes averaging over the homogeneous directions. Correlations 
among these velocity variation fields form the LES analogy of the Reynolds 
stresses or turbulent intensities. A characteristic turbulent intensity (vf) 1 ^ 2 is 
shown as a function of the normal coordinate x% in figure 21 in which we com- 
pare Leray and dynamic model predictions with filtered DNS data. We observe 
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Figure 4. Convergence and comparison of the profile of the streamwise turbulent intensity (ff) 1 ^ 2 
at t = 80 and A = 1/16. The filtered DNS data are marked with o and the convergence on 32 3 
(dash-dotted), 64 3 (dashed) and 96 3 (solid) is shown for the Leray model (top set of curves) and 
the dynamic model (lower set of curves) . 

a clear convergence of the LES predictions with increasing subgrid resolution 
A/h. The turbulent flow associated with either of these subgrid models ap- 
pears well resolved for A/h > 4 as observed above. The approximately grid- 
independent predictions arising from the Leray and dynamic models are seen to 
slightly underestimate the streamwise turbulent intensity in the developed tur- 
bulent regime. Compared with the dynamic model, the Leray predictions agree 
more closely with filtered DNS data. A similar level of agreement was found 
for the other turbulent intensities (v 2 ) 1 / 2 and the analogy of the Reynolds 
stress {viV2)- 



4.1.3 Resolved kinetic energy spectrum.. A more detailed assessment is 
obtained from the streamwise kinetic energy spectrum shown in figure 03 The 
dynamic model yields a significant under-prediction of the intermediate and 
smaller retained scales, particularly for the approximately grid-independent 
solution. The Leray predictions are much better. On coarse grids, the in- 
teraction with the spatial discretization method is seen to lead to a strong 
over-prediction of the smaller retained scales. However, at proper numerical 
subgrid resolution the situation improves considerably and the Leray model 
is seen to capture all scales with high accuracy. A slight under-prediction of 
the smaller scales and a small over-prediction of the large scales remains as 
systematic error in the grid-independent Leray solution. 
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Figure 5. Strcamwise kinetic energy spectrum E at t = 75: filtered DNS (o), Leray-modcl (32 3 : 
dash-dotted, 64 3 : solid, 96 3 : A), dynamic model (32 3 : dashed, 64 3 : dashed with o). A fixed 
filterwidth of £/16 was used and the Reynolds number Re = 50. 




Figure 6. Total (Tt), forward (1/), backward (T(,) energy transfer for the Leray model comparing 
three different grids. Thin lines correspond to 32 3 : Tf, dash-dotted, Tf dashed, Tt solid, thick lines 
correspond to 64 3 : T b dash-dotted, Tf dashed, T t solid and markers to 96 3 : T b o, Tf o, T t □. A 
fixed filterwidth of was used and Re = 50. 



4.1.4 Forward and backward scatter of energy.. The 'inner operations' 
of a subgrid model can be further characterised in terms of its dissipative or 
productive contributions to the evolution of the resolved kinetic energy E. 
This evolution is governed by 



(52) 
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where we introduced the viscous dissipation (T v ) and the total subgrid trans- 
fer (Tt). A positive total subgrid transfer T-t contributes to a reduction of the 
resolved kinetic energy and vice versa. The predicted kinetic energy evolu- 
tion, corresponding to a given subgrid model may be found by replacing the 
turbulent stress tensor by the assumed subgrid model. 

To further analyse the dynamical behaviour, one may split the total subgrid 
contribution Tj into a positive, i.e., forward scatter or dissipative, contribu- 
tion Tf and a negative, i.e., backward scatter or reactive contribution T . To 
formalise this splitting, we introduce 



Tf = - {uidjTij + \uidjTij\) dx. , T b = - {uidjTij - \uidjTij\) dx (53) 



In figure El we collected the forward, backward and total energy transfer con- 
tributions. The individual contributions grow in magnitude in the transitional 
and turbulent stages of the flow and are seen to be accurately captured nu- 
merically for A/h > 4. The Leray model is seen to contribute a significant 
backscatter component to the dynamics of the flow. Also, the total transfer 
becomes strictly positive beyond the transitional stages t > 50 which indicates 
that the Leray model on average adds to the dissipation of kinetic energy. 

Besides the performance of the Leray model at relatively low Reynolds num- 
ber, the predictions in the asymptotic high Reynolds number range are impor- 
tant in case one is interested in applying this subgrid model to flow problems 
of realistic complexity. In the next subsection we consider this in some more 
detail. 

4.2 High Reynolds number limit 

A particularly appealing property of Leray modelling is its robustness at very 
high Reynolds numbers. This is illustrated in figure where we collected spec- 
tra of the resolved kinetic energy in the turbulent regime at three different 
Reynolds numbers. The lower Re = 50 allows a corresponding direct numeri- 
cal simulation. However, at Re = 500 and Re = 5000 such a direct simulation 
is not feasible and flow predictions in this regime are possible only using LES. 
The robustness at very high Reynolds numbers of the Leray model is quite 
unique for a subgrid model without an explicit eddy-viscosity contribution. 
Although comparison with filtered DNS data is impossible here, we observe 
that the smoothed Leray dynamics is essentially captured as r = A/h > 4 |38| . 
The tail of the spectrum increases with Re, indicating a greater importance 
of small scale flow features. Improved subgrid resolution shows a reduction of 
these smallest scales, consistent with the reduced numerical error. At high Re 
the spectrum corresponding to the Leray model tends to contain a region with 
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approximately A;~ 5 / 3 behaviour, which is absent at Re = 50. 

The Leray prediction of the dependence of the streamwise velocity fluctu- 
ations on the Reynolds number is collected in figure Ufa). We observe that 
(vf) 1 / 2 decreases strongly with the Reynolds number, while the shape and 
width of the profiles remains quite similar for all Reynolds numbers. In this 
range, at fixed filter- width A, an increase in Re by a factor of 10 results in 
a similar reduction of the maximum of the streamwise velocity fluctuations. 
The subgrid transfer of resolved kinetic energy is illustrated in figure E^b). We 
observe that an increase in Re gives rise to an increase in the magnitude of 
each of the distinguished components to the evolution of E. The total transfer 
is positive in the developed stages, thereby contributing to an increased dissi- 
pation of energy as Re increases. The backscatter T appears to approach an 
asymptotic limit while the forward scatter Tf is still strongly increasing with 
Re. 

As was described in subsection 12. 31 the Leray model is part of the LANS— a 
regularisation. The LANS— a approach provides consistency of the large eddy 
simulation with the filtered Kelvin circulation theorem. In the next section 
we will investigate what consequences this important extension of the Leray 
model has for the LES capturing of turbulent mixing. 
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Figure 8. The Reynolds dependence of the streamwise velocity fluctuations are shown in (a): 
Re = 50 (dash-dotted), Re = 500 (dashed) and Re = 5000 (solid). In (b) the total (T t ), forward 
{Tt), backward (T b ) energy transfer for the Leray model are shown: thin lines correspond to 
Re = 50: Tj, dash-dotted, Tf dashed, Tt solid, thick lines correspond to Re = 500: TJ, dash-dotted, 
T f dashed, T t solid and markers to Re = 5000: T b o, T f o, T t □. A fixed filterwidth of £/16 and a 

grid with 96 3 cells was used. 

5 LANS— a improvements and limitations 

In this section we will illustrate some improvements over Leray and dynamic 
model predictions that arise when the LANS— a subgrid model is adopted. 
As explained in subsection 13.11 the LANS— a model appears to better rep- 
resent the small-scale variability that is contained in the turbulent flow and 
to correspond quite closely with filtered DNS data. This suggests also that 
derived macroscopic flow properties such as the resolved kinetic energy or the 
momentum thickness are better described using the LANS— a model. In sub- 
section 15.11 we will investigate a number of flow properties to quantify these 
improvements. A practical consequence of retaining more small-scale variabil- 
ity in the numerical solution is that the required resolution also needs to 
be increased. A discussion of these computational aspects and corresponding 
practical limitations for the LANS— a model is presented in subsection 15.21 

5.1 LANS— a prediction of turbulent mixing. 

In order to discuss the alterations in LES predictions that arise from adopting 
the LANS— a subgrid model, it is essential to distinguish between the findings 
associated with the approximately grid-independent solution and the conver- 
gence process toward this solution. This distinction should of course be made 
for every subgrid model but it is all the more relevant for the LANS— a pa- 
rameterisation since this subgrid model allows the existence of a significant 



Leray and LANS— a modelling of turbulent mixing 



35 



level of small-scale variability in the resolved flow. The LANS— a model corre- 
sponds very closely to filtered DNS snapshots, provided the subgrid resolution 
is adequate. Under-resolution of the LANS— a model will constitute a source 
for strong numerical contamination which is more characteristic of the spatial 
discretization that was used than a measure for the quality of the subgrid 
model. Therefore, in this subsection we will first turn to the prediction of the 
kinetic energy and momentum thickness obtained at high subgrid resolution 
and consider the convergence process afterwards. 




Figure 9. Resolved kinetic energy E (a) and momentum thickness <5 (b) at Re = 50 for A = £/16 

and three LES models: LANS— a (solid), Leray (dashed) and dynamic model (dash-dotted), 
compared with filtered DNS (solid circles) corresponding to the approximately grid-independent 

prediction at 96 3 . 



The evolution of the resolved kinetic energy and the momentum thickness 
is shown in figure In this figure the approximately grid-independent predic- 
tions obtained at a resolution 96 3 , i.e., r = 6 are collected. As was established 
in the previous section, we observe a strong improvement of the prediction 
of S using the Leray model compared to the dynamic model, whereas the re- 
solved kinetic energy is predicted at a similar level of accuracy. In addition, 
the LANS— a results are seen to constitute a further significant improvement 
for 5 and agree almost perfectly with the filtered DNS results. The resolved 
kinetic energy is also seen to be much better described by the LANS— a model; 
although the level of kinetic energy is slightly over-predicted, the decay rate 
dE/dt is very close to that seen in the filtered DNS data. For both E and 5 
the improvement is remarkable in the early laminar and transitional regime, 
which appears to set the stage for the accurate prediction of S in the developed 
turbulent flow. Only in the very late stages of the simulation is a deviation 
seen to develop in the prediction for S, and this may be due to a less than 
complete resolution of the smallest retained scales. These simulation results 
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clearly illustrate the improvements in the description of the physics of the flow 
by the LANS— a model which become available in case the spatial resolution 
is adequate. 
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Figure 10. Convergence and comparison of the profile of the streamwise turbulent intensity 
(v 2 ) 1 / 2 at t = 80 and with A = i/16. The filtered DNS data are marked with o. The convergence of 
the LANS— a predictions on 32 3 (dashed), 64 3 (dash-dotted) and 96 3 (solid) is shown in (a) and 
the approximately grid-independent results are shown in (b) for the Leray model (dashed), the 
dynamic model (dash-dotted) and the LANS— a model (solid). 



The prediction of the streamwise turbulent intensity is shown in figure I1U1 
We notice that under-resolution of the LANS— a model results in too high 
levels of turbulent intensity. Further visualisation of the LANS— a flow at a 
resolution of 32 3 showed snapshots which are dominated by an unphysically 
high level of small scale features. However, at appropriate spatial resolution we 
observe a fair approximation of a grid-independent solution that is seen to cor- 
respond very well with the filtered DNS data in figure lTOl' b). The LANS— a and 
Leray predictions are both accurate representations of the filtered DNS find- 
ings with a slight improvement of the predicted maximal turbulent intensities 
near the centerline due to the LANS— a model. There are certainly consider- 
able improvements over the results that are obtained when using the dynamic 
model. To put this further into perspective, an earlier study by Vreman et 
al. established that the dynamic model was among the more accurate 
models compared to the Bardina similarity model, the nonlinear or gradient 
model and other dynamic (mixed) models. 

The type of contribution of the LANS— a model to the evolution of the 
resolved kinetic energy is displayed in figure 111! Although the convergence 
toward the grid-independent solution is not as dramatic as was observed in 
relation to the Leray model in figure the strong backscatter effects of the 
LANS— a model are well established. We observe that the total energy transfer 
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Figure 11. Total (Ti), forward (TV), backward (TJ,) energy transfer for the LANS— a model 
comparing three different grids in (a): thin lines correspond to 32 3 : T b dash-dotted, TV dashed, Tt 
solid, thick lines correspond to 64 3 : T b dash-dotted, TV dashed, Tt solid and markers to 96 3 : T;, o, 
TV. o, T t □. In (b) we compare the fine-grid results for the Leray and LANS— a model in (b): Leray 
(thick lines) and LANS— a (thin lines): T t (solid), Tf (dashed) and T b (dash-dotted). A fixed 
filterwidth of £/16 was used and the Reynolds number Re = 50. 

T t is negative which is an essential difference compared to the Leray model. 
This is best illustrated in figure UTT b) where the kinetic energy transfers of the 
Leray and LANS— a models are compared. We observe that the levels of for- 
ward and backward scatter are considerably more prominent in the LANS— a 
model. Moreover, the Leray model displays a slight negative total energy trans- 
fer only in the earlier stages of the flow while the LANS— a model shows a 
total backscatter during the entire evolution. This is the underlying reason for 
the earlier observed higher level of variability in the LANS— a flow. 

The increased small-scale contributions in the LANS— a predictions require 
appropriate spatial resolution. The corresponding resolution requirements are 
discussed in some more detail next. 



5.2 Resolution requirement for the LANS— a model. 

To more clearly express the deviations in the momentum thickness compared 
to filtered DNS data, one may choose to amplify the total simulation errors 
that arise. We observe that the mixing layer displays a nearly linear, self- 
similar dependence of 8 as function of time. To test this linearity and hence 
the convergence process more precisely, we introduce 
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Figure 12. Convergence of momentum thickness <5 at Re = 50 for A = £/16 and three LES models: 
LANS— a (solid), Leray (dashed) and dynamic model (dash-dotted), compared with filtered DNS 
(solid circles). We plotted e(t) = & LE s(t) - %jvs (**)(*/**) for ** = 100: 323; no markers, 64 3 : 

open circles and 96 3 : thick lines. 



and use as reference time t* = 100 which corresponds to the final time used in 
the simulation. In figure fH21 we consider the convergence of the predictions as 
a function of spatial resolution for the three subgrid models adopted in this 
paper. We observe that the dynamic model results are less accurate; but the 
solution is numerically captured to its full potential already at 32 3 . The Leray 
predictions are seen to require 64 3 cells in order to attain their full potential 
for accurate prediction. The results of the LANS— a model are slightly more 
sensitive and approach grid-independence only at 96 3 . This test of linearity of 
5 also illustrates the consequences of under-resolution. Although the LANS— a 
predictions are the most accurate among the models at proper subgrid reso- 
lution, the effects of numerical contamination at insufficient resolution can be 
strong enough to lose most of this potential, e.g., on a grid with 32 3 cells. This 
pattern of dependency of its predictions on the spatial resolution when the 
LANS— a model is under-resolved was also observed for other flow properties, 
such as the decay of the resolved kinetic energy. 

From these numerical illustrations we may infer a rule for the spatial res- 
olution requirements of the LANS— a model. As discussed in subsection 12.31 
the LANS— a model derives its name from the length-scale parameter a. By 
considering the Taylor expansion of the top-hat filter-operation and compar- 
ing this with the inverse Helmholtz operator one arrives at the identification 
a ~ A/5 30 . Combined with the observation made above that a reliable treat- 
ment of the LANS— a model requires a subgrid resolution r 6 leads to the 
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conclusion that LES of turbulent mixing based on the LANS— a model requires 
a « h. 

6 Concluding remarks 
6.1 Summary of results 

In this paper, we proposed to consider the mathematical approaches for reg- 
ularisation of the Navier-Stokes equations as models for parameterising the 
effects of the unresolved length-scales in large eddy simulation of turbulent 
mixing flow. Two related regularisation principles were considered, namely, the 
Leray approach and the LANS— a approach. These regularisation approaches 
produced accurate LES predictions that do not depend on additional ad hoc 
implementation steps, such as those required in implementing the dynamic 
procedure. The LANS— a model is a deformation of the Leray model which 
recovers the Kelvin circulation theorem of the Navier-Stokes equations, but for 
a material loop moving with a filtered transport velocity. In a simple Fourier 
analysis, these regularisation models were compared to Bardina's similarity 
model and to the nonlinear tensor-diffusivity model. Expansions of the turbu- 
lent stress tensor provided a point of reference for these comparisons. A more 
meaningful assessment of the quality of these subgrid models was obtained by 
comparing them in numerical simulations of turbulent transport in a temporal 
mixing layer. At the level of instantaneous solutions, the regularisation mod- 
els were shown to correspond much better with filtered DNS data than has 
been seen for other, more traditional subgrid models in literature. This general 
impression also translates into a better prediction of various flow properties, 
ranging from mean flow to spectral quantities. 

The Leray model was shown to provide a representation of turbulent mixing 
that is more accurate in many respects than the predictions obtained using the 
dynamic eddy-viscosity model. Moreover, on any simulation grid the computa- 
tional effort required for the Leray model is considerably lower than required 
for the dynamic model. The resolved kinetic energy was found to be slightly 
overestimated by the Leray model, but was otherwise very similar to that for 
the dynamic model. Important improvements were observed in the capturing 
of the momentum thickness and the velocity fluctuation profiles. Particularly, 
the intermediate and smaller resolved scales in the turbulent regime were much 
better represented by the Leray model, compared to the dynamic model. In 
addition, the intricacies of turbulent energy dynamics are contained more fully 
in the Leray model, because it allows both forward and backward scatter of 
energy. Finally, the Leray model was found to be robust at very high Reynolds 
number and the prediction of the dominant part of the resolved kinetic energy 
spectrum was found to approach a —5/3 behaviour with increasing Reynolds 
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number. 

Compared to the Leray and the dynamic subgrid models, the LANS— a 
model was shown to yield significant improvements, provided the spatial reso- 
lution is such that a ~ h ~ A/5. The improvements associated with the grid- 
independent solution displayed a rather strong deterioration in cases where 
the resolution was not adequate. It was observed that at subgrid resolution 
ratio r = A/h = 2 the accuracy of the LANS— a predictions decreased to 
being the same as that of the Leray and dynamic models. The requirement of 
adequate subgrid resolution posed some limitations on the practical use of the 
regularisation models. We concluded for the Leray model that r > 4 and for 
the LANS— a model that r > 6 constitute reliable values for the specification 
of the spatial resolution. 

6.2 Comparison of the regularisation models and outlook 

The Leray model displays excellent robustness with increasing Reynolds num- 
ber. This feature allows one to apply the Leray model accurately at reasonable 
computational costs and under flow-conditions that are well outside current 
DNS capabilities. However, the LANS— a model yields solutions with more 
realistic variability, corresponding better to the filtered DNS results than for 
the Leray model. Thus, a trade-off emerges between these two models. The 
solutions of the LANS— a model may more accurately represent the effects 
of intermittency in turbulence than the less-variable solutions of the Leray 
model. However, the LANS— a model is less robust and its application to flow 
at high Reynolds numbers is not as straightforward as with the Leray model. 
Further investigation of this trade-off may lead to interesting developments in 
the comparison of the time-dependent solutions of these two models. 

A convenient benefit of the regularisation approach to turbulence modelling 
is that it enables one to derive the implied small-scale treatment from the 
underlying regularisation principle. This yields a systematic closure of the 
equations whose analysis allows an extension in which the filter width A is 
determined dynamically by the evolving flow. The evolving filter-width may 
even be anisotropic [401141] . The application of this self-adaptive modelling 
approach in a spatially developing mixing layer and, more importantly, in 
near wall turbulence is a topic of current research. 

Acknowledgement 

BJG gratefully acknowledges support from the Turbulence Working Group 
(TWG)at the Centre for Non-Linear Studies (CNLS) at Los Alamos National 
Laboratory which facilitated an extended research visit in 2004 and allowed 



Leray and LANS— a modelling of turbulent mixing 



41 



many fruitful discussions with members of TWG. Simulations were performed 
at SARA and made possible through grant SC-244 of the Dutch National 
Computing Foundation (NCF). 



References 

Ghosal, S.: 1999. Mathematical and physical constraints on large eddy simulation of turbulence. 
AIAA J., 37, 425. 

Spcziale, C.G.: 1991. Analytical Methods for the Development of Reynolds-Stress Closures in 
Turbulence, Ann. Rev. Fl. Mech. 23, 107. 

Horiuti, K.: 2003. Roles of non-aligned eigenvectors of strain-rate and subgrid-scale stress tensors 
in turbulence generation, J. Fluid Mech. 491, 65. 

Vreman, A.W., Geurts, B.J., Kuerten, J.G.M.: 1994. Realisability conditions for the turbulent 
stress tensor in large eddy simulation. J. Fluid Mech., 278, 351. 
Germano, M.: 1992. Turbulence: the filtering approach. J. Fluid Mech., 238, 325. 
Geurts, B.J.: 1997. Inverse modelling for large eddy simulation. Phys. of Fluids, 9, 3585. 
Kuerten, J.G.M., Geurts, B.J., Vreman, A.W., Germano, M.: 1999. Dynamic inverse modelling 
and its testing in large eddy simulations of the mixing layer. Phys. of Fluids, 11, 3778. 
Sagaut, P.: 2001. Large eddy simulation for incompressible flows; an introduction. Scientific 
Computation. Springer Verlag. 

Smagorinsky, J.: 1963. General circulation experiments with the primitive equations. Mon. 
Weather Rev., 91, 99. 

Berselli, L.C., Galdi, G.P., Ilicscu, T., Layton, W.J.: 2002. Mathematical Analysis for the Ra- 
tional Large Eddy Simulation model, Math. Mod. Meth.Appl.Sci., 12, 1131. 
Vreman, A.W., Geurts, B.J., Kuerten, J.G.M.: 1997. large eddy simulation of the turbulent 
mixing layer. J. Fluid Mech., 339, 357. 

Leray, J.: 1934. Sur les movements d'un fluide visqueux remplaissant l'espace. Acta Mathematica, 
63, 193. 

Foias, C, Holm, D.D., Titi, E.S.: 2001. The Navier-Stokes-alpha model of fluid turbulence. 
Physica D, 152, 505. 

Cheskidov, A., Holm, D.D., Olson, E.J., Titi, E.S.: 2005. On a Leray-alpha model of turbulence. 
Phys. Fluids To appear. 

Holm, D.D., Marsden, J.E., Ratiu, T.S.: 1998 The Euler-Poincare equations and semidirect 
products with applications to continuum theories. Adv. in Math., 137, 1. 

Germano, M., Piomelli U., Moin P., Cabot W.H.: 1991. A dynamic subgrid-scale eddy viscosity 
model. Phys. of Fluids, 3, 1760 

Bos, F. van der, Geurts, B.J.: 2005. Commutator errors in the filtering approach to large eddy 
simulation. Phys. of Fluids, 17, 035108 

Geurts, B.J.: 2003. Elements of direct and large eddy simulation. Edwards Publishing, Inc. 
Meneveau, C, Katz, J.: 2000. Scalc-invariance and turbulence models for large eddy simulation. 
Annu. Rev. Fluid Mech., 32, 1. 

Stolz, S., Adams, N.A.: 1999. An approximate deconvolution procedure for large eddy simulation. 
Phys. of Fluids, 11, 1699. 

Kolmogorov, A.N.: 1991. The local structure of turbulence in an incompressible fluid with very 
large Reynolds numbers. Proc. Roy. Soc. London, 434, 1890 (First published: 1941, Dokl. Akad. 
Nauk SSSR, 30, 301). 

Frisch, U.: 1995. Turbulence: the legacy of A.N. Kolmogorov. Cambridge University Press, Cam- 
bridge. 

Bardina, J., Ferziger, J.H., Reynolds, W.C.: 1984. Improved turbulence models based on LES of 
homogeneous incompressible turbulent flows. Dept. Mech. Engg. Report No. TF-19, Stanford. 
Vreman, A.W., Geurts, B.J., Kuerten, J.G.M.: 1994. On the formulation of the dynamic mixed 
subgrid-scale model. Phys. of Fluids, 6, 4057. 

Meneveau, C, Lund, T.S., Cabot, W.H.: 1996. A Lagrangian dynamic subgrid-scale model of 
turbulence. J. FLuid Mech., 319, 353. 

Lilly, D.K.: 1992. A proposed modification of the Germano subgrid-scale closure method. Phys. 
of Fluids A, 4, 633. 



42 Leray and LANS— a modelling of turbulent mixing 

[27] Vreman, A.W.: 1995. Direct and large eddy simulation of the compressible turbulent mixing 

layer. Ph.D. Thesis, University of Twente. 
[28] Gcurts, B.J., Holm, D.D.: 2003. Rcgularisation modelling for large eddy simulation. Phys. of 

Fluids, 15, L13. 

[29] Holm, D.D.: 1999. Fluctuation effects on 3D Lagrangian mean and Eulerian mean fluid motion. 
Physica D, 133, 215. 

[30] Gcurts, B.J., Holm, D.D.: 2002. Alpha-modelling strategy for LES of turbulent mixing. Turbulent 
Flow Computation. Eds: D. Drikakis, B.J. Geurts. Fluid mechanics and its applications 66, 
Kluwcr Academic Publishers, 237. 

[31] Domaradzki, J. A., Holm, D.D.: 2001. Navier-Stokes-alpha model: LES equations with nonlin- 
ear dispersion. Modern simulation strategies for turbulent flow. Edwards Publishing, Ed. B.J. 
Geurts. 107. 

[32] Gcrmano, M.: 1986. Differential filters for the large eddy numerical simulation of turbulent flows. 
Phys. of Fluids, 29, 1755. 

[33] Winckelmans, G.S., Jeanmart, H., Wray, A. A., Carati, D., Geurts, B.J.: 2001. Tcnsor-diffusivity 
mixed model: balancing reconstruction and truncation. Modern simulation strategies for turbu- 
lent flow. Edwards Publishing, Ed. B.J. Geurts. 85. 

[34] Vreman, A.W., Gcurts, B.J., Kuerten, J.G.M.: 1996. Large eddy simulation of the temporal 
mixing layer using the Clark model. TCFD, 8, 309. 

[35] Cao, C, Holm, D., Titi, E.S.: 2004: On the Clark— a model of turbulence: global regularity and 
long-time dynamics, Journal of Turbulence, Submitted. 

[36] Jameson, A.: 1983. Transonic flow calculations. MAE-Report 1651, Princeton University. 

[37] Geurts, B.J., Kuerten, J.G.M.: 1993. Numerical aspects of a block-structured flow solver. J. Eng. 
Math., 27, 293. 

[38] Geurts, B.J., Froehlich, J.: 2002. A framework for predicting accuracy limitations in large eddy 
simulation. Phys. of fluids, 14, L41. 

[39] Meyers, J., Gcurts, B.J., Baelmans, M.: 2003. Database analysis of errors in large eddy simula- 
tion. Phys. of Fluids, 15, 2740. 

[40] Holm, D.D.: 1999 Fluctuation effects on 3D Lagrangian mean and Eulerian mean fluid motion. 
Physica D, 133, 215. 

[41] Holm, D.D.: 2002 Lagrangian averages, averaged Lagrangians, and the mean effects of fluctua- 
tions in fluid dynamics." Chaos 12, 518. 



